Loading up the data for some preliminary analyses of the binary climb/no-climb categories.
You can see the indices I add in the code below. Note that the indices are all caps, and the linear measurements have lowercase letters after the first.
dat <- read_sheet("https://docs.google.com/spreadsheets/d/1-eknhyZ1JNnXqhg2kViyzVntC8NGZvILQX-aQQb1Jvk/edit#gid=325036460", na = c("NA", "?", "")) %>%
select(!NOTES) %>%
# Recode Ordinal Rankings
mutate(Loc_Ord = case_when(Loc_mode_Ordinal == "G" ~ 1,
Loc_mode_Ordinal == "A" ~ 2,
Loc_mode_Ordinal == "Sc" ~ 3,
Loc_mode_Ordinal == "T" ~ 4,
Loc_mode_Ordinal == "Is" ~ 5,
Loc_mode_Ordinal == "Sf" ~ 5,
Loc_mode_Ordinal == "Ss" ~ 6,
TRUE ~ NA),
Loc_Ord = as.ordered(Loc_Ord),
Loc_bin = case_when(Loc_mode_Bindary == "Ground" ~ 0,
Loc_mode_Bindary == "Tree" ~ 1,
TRUE ~ NA
),
# Loc_bin = as.factor(Loc_bin),
Loc_mode_Categorical = as.factor(Loc_mode_Categorical),
log_Mass = log(Mass_grams)) %>%
relocate(Loc_bin, .after = Loc_mode_Bindary) %>%
relocate(Loc_Ord, .after = Loc_mode_Ordinal) %>%
relocate(log_Mass, .before = Skl) %>%
#################
#Calculate Indices!
#################
mutate(SI = Sh / Sl, # Scapular Index
HRI = Hsw / Hl, # Humeral Robustness Index
HPI = Hpw / Hl, # Humeral Proximal Index
HEB = Hdw / Hl, # Humeral Epicondyle Breadth
HHRI = Hhl / Hl, # Humeral Head Robustness Index
HHW = Hhw / Hhl, # Humeral Head Shape Index
DI = Hdcw / Hsw, # Deltopectoral Crest Index
OLI = Uol / Ul, # Olecranon Process Length Index
BI = Rl / Hl, # Brachial Index
IM = (Hl+Ul)/(Fl+Tl), # Intermembral Index
PRTI = Mcl/(Hl+Rl), # Palm Robustness Index
MRI = Mcw / Mcl, # Metacarpal Robustness
MANUS = Ppl / Mcl, # MANUS index
MANUS2 = (Ppl+Ipl)/Mcl, # MANUS index with intermed. phalanx
IRI = Fgh / Fl, # Gluteal Index
FRI = Fsw / Fl, # Femoral Robustness
FEB = Fdw / Fl, # Femoral Epicondyle Breadth
CI = Tl / Fl, # Crural Index
TRI = Tmw / Tl, # Tibial Robustness Index
#ANR = Anl / Al, # Astragular Neck Robustness Index
#CAR = Cal / Cl, # Calcaneal Robustness Index
IRI = Il / Pel, # Illium Robustness Index
PR = Il / Isl, # Pelvic Index
PES = Pppl / Mtl, # PES INdex
PES2 = (Pppl+Pipl)/Mtl # PES with intermediate Phalanx
) %>%
mutate_at(vars(17:71), log) %>%
mutate_at(vars(16:93), scale2)
What does the missing data look like?
You can scroll through the table below to see
| measurement | count_missing | percent_missing |
|---|---|---|
| Pdpw | 324 | 0.773 |
| Pipw | 278 | 0.663 |
| Pdpl | 277 | 0.661 |
| Dpw | 263 | 0.628 |
| Jl | 249 | 0.594 |
| Anl | 239 | 0.570 |
| Atw | 239 | 0.570 |
| Cal | 239 | 0.570 |
| Ccw | 238 | 0.568 |
| Csw | 238 | 0.568 |
| Ctl | 238 | 0.568 |
| Ctw | 238 | 0.568 |
| Al | 237 | 0.566 |
| Pppw | 230 | 0.549 |
| Skl | 228 | 0.544 |
| Mtw | 223 | 0.532 |
| Ipw | 217 | 0.518 |
| Dpl | 216 | 0.516 |
| PES2 | 194 | 0.463 |
| Pipl | 193 | 0.461 |
| Fbdw | 189 | 0.451 |
| Fbmw | 187 | 0.446 |
| Fgh | 187 | 0.446 |
| Fhd | 187 | 0.446 |
| HHRI | 185 | 0.442 |
| HHW | 185 | 0.442 |
| Hhl | 185 | 0.442 |
| Hhw | 185 | 0.442 |
| Ppw | 158 | 0.377 |
| Cl | 154 | 0.368 |
| Fbpw | 153 | 0.365 |
| MRI | 152 | 0.363 |
| Mcw | 152 | 0.363 |
| DI | 144 | 0.344 |
| Hdcw | 144 | 0.344 |
| Tdw | 141 | 0.337 |
| Tpw | 138 | 0.329 |
| Fbl | 137 | 0.327 |
| IRI | 137 | 0.327 |
| Isl | 137 | 0.327 |
| PR | 137 | 0.327 |
| Pel | 137 | 0.327 |
| Il | 136 | 0.325 |
| HPI | 135 | 0.322 |
| Hpw | 135 | 0.322 |
| Ipl | 117 | 0.279 |
| MANUS2 | 117 | 0.279 |
| PES | 95 | 0.227 |
| Pppl | 94 | 0.224 |
| Mtl | 88 | 0.210 |
| MANUS | 23 | 0.055 |
| Ppl | 23 | 0.055 |
| Mcl | 17 | 0.041 |
| PRTI | 17 | 0.041 |
| log_Mass | 10 | 0.024 |
| TRI | 3 | 0.007 |
| Tmw | 3 | 0.007 |
| CI | 2 | 0.005 |
| FEB | 2 | 0.005 |
| Fdw | 2 | 0.005 |
| IM | 2 | 0.005 |
| Tl | 2 | 0.005 |
| FRI | 1 | 0.002 |
| Fl | 1 | 0.002 |
| Fsw | 1 | 0.002 |
| SI | 1 | 0.002 |
| Sh | 1 | 0.002 |
| Sl | 1 | 0.002 |
| BI | 0 | 0.000 |
| HEB | 0 | 0.000 |
| HRI | 0 | 0.000 |
| Hdw | 0 | 0.000 |
| Hl | 0 | 0.000 |
| Hsw | 0 | 0.000 |
| OLI | 0 | 0.000 |
| Rl | 0 | 0.000 |
| Ul | 0 | 0.000 |
| Uol | 0 | 0.000 |
Preliminary data analysis, looping over all of the variables to see which ones do a good job predicting the binary tree vs. no-tree categorization
These are very preliminary, and the results will become more nuanced and probably more confusing, but it can at least give us a sense of which measurements are informative of climbing.
Here are all the model plots. On the y axis, 1 is TREE, 0 is NO TREE. The x axis is the phenotype, mean centered on 0 and scaled to a standard deviation of 1. All of the linear measurements are log transformed prior to mean-centering. All the models include log_mass as a variable, meaning that they are “size corrected” representations of the effect of the phenotype on climbing. What we are looking for is a slope that ranges across the whole y axis, meaning that it touches the 1 and 0, and has a steep slope (either up or down). To interpret, look at the 3rd plot, Hl. As Hl increases, the probability of being TREE increases, or, the longer the humerus, the more likely to climb!
Here are some standouts: (remember, these are log-transformed and size-corrected effect sizes)
for(i in fit_list2){
plots <- plot(conditional_effects(i), plot=F, points = T)[[1]]
print(plots + theme_bw())
}
Repeating above, but with Species as a group-level effect. There are a total of 243 species, and 57 have more than one sample.
Plots
for(i in fit_list3){
plots <- plot(conditional_effects(i), plot=F, points = T)[[1]]
print(plots + theme_bw())
}